/*
 * ExtendedTreeSolver.hpp
 *
 * Copyright 2012 fabian
 *
 * This file is part of Recognition-Procedures-for-Boolean-Functions.
 *
 * Recognition-Procedures-for-Boolean-Functions is free software: you
 * can redistribute it and/or modify it under the terms of the
 * GNU General Public License as published by the Free Software Foundation,
 * either version 3 of the License, or (at your option) any later version.
 *
 * Recognition-Procedures-for-Boolean-Functions is distributed in the hope
 * that it will be useful, but WITHOUT ANY WARRANTY; without even the implied
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General
 * Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with Recognition-Procedures-for-Boolean-Functions.
 *
 * If not, see <http://www.gnu.org/licenses/>.
 */

#ifndef EXTENDEDTREESOLVER_HPP_
#define EXTENDEDTREESOLVER_HPP_

#include <string>
#include <iostream>
#include <utility>
#include <algorithm>

#include "bfr/LPBConversionException.hpp"
#include "bfr/lpbCombinatorial/SplittingTree.hpp"
#include "bfr/lpbCombinatorial/SolverInformation.hpp"
#include "bfr/lpbCombinatorial/CoefficientSelectionException.hpp"
#include "bfr/LPB.hpp"

namespace bfr {

template<class T>
std::ostream& operator<<(std::ostream &out, const std::pair<T, T> &p) {
        out << "(" << p.first << ", " << p.second << ")";
        return out;
}

union RepairInputContent {
        SplittingTreeNode *node;
        // 0 -> false node, 1 -> true node
        char trueType;
};

struct RepairInput {
    public:
        explicit RepairInput(bool is_final)
        : final(is_final) { }
        bool isFinal() const {
            return final;
        }
        RepairInputContent getContent() const {
            return content;
        }

        void setContent(const RepairInputContent &c) {
            content = c;
        }

        void setFinal(bool f) {
            final = f;
        }

    private:
        bool final;
        RepairInputContent content;
};

struct RepairResult {
    public:
        explicit RepairResult(bool s=false, bool il=false,
                lpb_coefficient_t b=-1, size_t k_val=-1)
        : success(s), isLower(il), bound(b), k(k_val) { }

        bool success;
        bool isLower;
        lpb_coefficient_t bound;
        size_t k;
};

template<class ChooserT>
class ExtendedTreeSolver {
    public:
        explicit ExtendedTreeSolver(const ChooserT &c=ChooserT())
        : chooser(c), information() { }

        ~ExtendedTreeSolver() { }

        std::pair<RepairInput, RepairInput> createResults(const RepairInput &in) {
            RepairInput first(true);
            RepairInput second(true);
            RepairInputContent c1;
            RepairInputContent c2;
            SplittingTreeNode *n = in.getContent().node;
            if (in.isFinal()) {
                if (in.getContent().trueType == 0) {
                    c1.trueType = 0;
                    c2.trueType = 0;
                }
                else {
                    c1.trueType = 1;
                    c2.trueType = 1;
                }
            } else if (n->getUpperChild() == 0) { // false node with cutTree
                c1.trueType = 0;
                c2.trueType = 0;
            } else if (n->getLowerChild() == 0) { // true node with cutTree
                c1.trueType = 1;
                c2.trueType = 1;
            } else {
                if ((n->getUpperChild()->getUpperChild() == 0)) {
                    if (n->getUpperChild()->getDNF()->size() == 0) {
                        c1.trueType = 0;
                    } else {
                        assert(n->getUpperChild()->getDNF()->size() == 1);
                        c1.trueType = 1;
                    }
                } else {
                    first.setFinal(false);
                    c1.node = n->getUpperChild();
                }
                
                if (n->getLowerChild()->getLowerChild() == 0) {
                    if (n->getLowerChild()->getDNF()->size() == 0) {
                        c2.trueType = 0;
                    } else {
                        assert(n->getLowerChild()->getDNF()->size() == 1);
                        c2.trueType = 1;
                    }
                } else {
                    second.setFinal(false);
                    c2.node = n->getLowerChild();
                }
            }
            first.setContent(c1);
            second.setContent(c2);
            return std::make_pair(first, second);
        }

        /**
         * Choose the degree in the interval \f$(a, b]\f$.
         *
         * Tries to choose the degree of the LPB in the interval
         * \f$(a, b]\f$. If this is not possible an exception is thrown.
         *
         * BIG TODO
         * We just test if there is an integer solution, but because a and
         * b are integers as well this should be ok.
         *
         * Can someone think about it? ;)
         *
         * @param a First value of the interval.
         * @param b Second value of the interval.
         */
        lpb_coefficient_t chooseDegree(lpb_coefficient_t a,
                lpb_coefficient_t b)
        const throw(LPBConversionException) {
            if ((a == coeff_pos_infinity) || (b == coeff_neg_infinity))
                throw LPBConversionException(errorMessage(a, b));
            else if (a == coeff_neg_infinity) {
                if (b == coeff_neg_infinity) {
                    throw LPBConversionException(errorMessage(a, b));
                }
                else if (b == coeff_pos_infinity) {
                    return 0;
                }
                else {
                    // b is a number
                    return b;
                }
            }
            else {
                // a is a number
                if (b == coeff_pos_infinity) {
                    return a + 1;
                }
                else {
                    // a and b are numbers
                    lpb_coefficient_t smallestSol = a + 1;
                    if (smallestSol > b) {
                        // TODO: what if we can choose the coefficient,
                        // but not as an integer?
                        if (a < b) {
                            std::cerr << "BUG: See TreeToLPB chooseDegree."
                                    << "Degree just can't be an integer."
                                    << std::endl;
                        }
                        throw LPBConversionException(errorMessage(a, b));
                    }
                    else {
                        return smallestSol;
                    }
                }
            }
        }

        void setRepairFlag(size_t flag) {
            information.setRepairFlag(flag);
        }
        size_t getRepairFlag() {
            return information.getRepairFlag();
        } 

        LPB *solve(SplittingTree *tree) throw(LPBConversionException) {
            // actually we do not start with column index tree.size() - 1
            // because then we must use while (k >= 0) what is not good
            // for size_t because it will always be >= 0 
	        size_t termination_hack = 0;
            size_t k = tree->size();
            //std::cout << "tree size " << k << std::endl;
            //std::cout << "Number of final true nodes = " << tree->numFinalTrueNodes() << std::endl;
            //std::cout << "Number of clauses = " << tree->getRoot()->getDNF()->size() << std::endl;
            //std::cout << "Symmetry saving: final true nodes / number of clauses = " << std::endl;  
	    //std::cout << (float)tree->numFinalTrueNodes() / (float)tree->getRoot()->getDNF()->size() << std::endl;
            const size_t numVariables = tree->size() - 1;
            information.init(*tree);
            while (k > 0) {
	      //std::cout << "k = " << k << std::endl;
                information.resetIntervalFactor(k - 1);
		//std::cout << "Before handleColumn " << std::endl;
                interval_t res = tree->handleColumn(k - 1, chooser, information);
		//std::cout << "After handleColumn " << std::endl;
                const std::pair<lpb_coefficient_t, lpb_coefficient_t> bounds = information.getRepairBound(k - 1);
		bool lBoundFromRepair = false;
		bool uBoundFromRepair = false;
                lpb_coefficient_t max = std::max(res.first, bounds.first);
		if (max > res.first) { lBoundFromRepair = true; }//JGS: provisional!!
		//std::cout << "max(" << res.first << ", " << bounds.first << ") = " << max << std::endl;
                lpb_coefficient_t min = std::min(res.second, bounds.second);
		if (min < res.second) { uBoundFromRepair = true; }//JGS: provisional!!
		//std::cout << "min(" << res.second << ", " << bounds.second << ") = " << min << std::endl;
                if (k != tree->size()) {
                    const lpb_coefficient_t lastCoeff = information.getCoefficient(k);
                    max = std::max(max, lastCoeff - 1);
                }
                if ((max + 1) == min) {
                    // yes, we have a conflict
                    // use the information to solve the conflict
                    //std::cout << "doubling the coefficients etc." << std::endl;
                    interval_t newValues =
                            information.solveConflict(k - 1, max, min);
                    max = newValues.first;
                    min = newValues.second;
                }
                try {
                    // different from java version, but should be correct
                    if (k != 1) {
                        // TODO k - 1 correct?
                        const lpb_coefficient_t coeff =
                                chooser.chooseCoefficient(max, min,
                                        tree->getContext(),
                                        information, k - 1);
                        //std::cout << "Choosing between " << max << ", " << min << std::endl;
                        //std::cout << "set " << k << " to " << coeff << std::endl;
                        information.setCoefficient(k - 1, coeff);
                    }
                } catch (const CoefficientSelectionException &e) {
                    // throw LPBConversionException(e.what());
            //BEGIN mic mac stuff
                    SplittingTreeNode *minCulprit = tree->getMinCulprit();
                    SplittingTreeNode *maxCulprit = tree->getMaxCulprit();

                    RepairInput mic(true);
                    RepairInput mac(true);
                    RepairInput micSibling(true);
                    RepairInput macSibling(true);

                    //JGS: moved these 4 lines up
                    //the next 20 something lines should be brushed up for readability
                    RepairInputContent c1;
                    RepairInputContent c2;
                    RepairInputContent c3;
                    RepairInputContent c4;

                    //JGS: the next 4 elses added by me
                    if (minCulprit->getUpperChild() != 0) {
                        mic.setFinal(false);
                        c1.node = minCulprit;
                    } else {
                        if (minCulprit->getDNF()->size() == 0) {
                            c1.trueType = 0; } else {
                                c1.trueType = 1;}
                    }

                    if (maxCulprit->getUpperChild() != 0) {
                        mac.setFinal(false);
                        c3.node = maxCulprit;
                    } else {
                        if (maxCulprit->getDNF()->size() == 0) {
                            c3.trueType = 0; } else {
                                c3.trueType = 1;}
                    }
                    
                    
                    if (minCulprit->getLowerParent()->getLowerChild()->getUpperChild() != 0) {
                        micSibling.setFinal(false);
                        c2.node = minCulprit->getLowerParent()->getLowerChild();
                    } else {
                        if (minCulprit->getLowerParent()->getLowerChild()->getDNF()->size() == 0) {
                            c2.trueType = 0; } else {
                                c2.trueType = 1;}
                    }

                    if (maxCulprit->getLowerParent()->getLowerChild()->getUpperChild() != 0) {
                        macSibling.setFinal(false);
                        c4.node = maxCulprit->getLowerParent()->getLowerChild();
                    } else {
                        if (maxCulprit->getLowerParent()->getLowerChild()->getDNF()->size() == 0) {
                            c4.trueType = 0; } else {
                                c4.trueType = 1;}
                    }
                    
                    mic.setContent(c1);
                    micSibling.setContent(c2);
                    mac.setContent(c3);
                    macSibling.setContent(c4);
		    //END mic mac stuff

                    //~ std::cout << "input final?" << std::endl;
                    //~ std::cout << mac.isFinal() << ", " << macSibling.isFinal() << ", "
						//~ << mic.isFinal() << ", " << micSibling.isFinal() << std::endl;
                    RepairResult r;
                    if (uBoundFromRepair) {
		      //std::cout << "!1" << std::endl;
                        if (lBoundFromRepair) {
                            //std::cout << "Both bounds are from repair. No idea what to do." << std::endl;
                        }
                        //std::cout << "uBoundFromRepair: call monorepair." << std::endl;
                        termination_hack++;//JGS: a barbaric way of ensuring termination of monorepair
                        if (termination_hack <= 100) {
                            information.setRepairFlag(1);
                            r = monorepair(mac, macSibling, k - 1, numVariables, min, false);
                        }
                    } else {
                        if (lBoundFromRepair) {
                            //std::cout << "lBoundFromRepair: call monorepair." << std::endl;
                            termination_hack++;//JGS: a barbaric way of ensuring termination of monorepair
                            if (termination_hack <= 100) {
                                information.setRepairFlag(1);
                                r = monorepair(mic, micSibling, k - 1, numVariables, max, true);
                            }
                        }
                        else {
                            information.setRepairFlag(1);
                            //std::cout << "call repair" << std::endl;
                            r = repair(mac, macSibling, mic, micSibling, k - 1, numVariables);
                        }
                    }
                    size_t shift_repair_col = 0;
                    if (!r.success) {
                        throw LPBConversionException(e.what());
                    }
                    else {
                        if (r.isLower) {
                            if (information.setLowerRepairBound(r.k, r.bound)) {
                                //std::cout << "setting lower bound in column " << r.k+1 << " to " << r.bound << std::endl;
                            }
                            else {
                                throw LPBConversionException("Lower bound could not be increased.");
                            }
                        }
                        else {
                            //if the upper bound computed by repair is smaller than the coefficient to the right of it,
                            //then the bound should instead be imposed on that coefficient
                            //std::cout << "r.bound = " << r.bound << std::endl;
                            //std::cout << "Coefficient = " << information.getCoefficient(r.k + shift_repair_col + 1) << std::endl;
                            while (r.bound <= information.getCoefficient(r.k + shift_repair_col + 1)) {
                                //std::cout << "Coefficient = " << information.getCoefficient(r.k + shift_repair_col + 1) << std::endl;
                                ++shift_repair_col; 
                                //std::cout << "shift_repair_col = " << shift_repair_col << std::endl; 
                            }
                            if (information.setUpperRepairBound(r.k + shift_repair_col, r.bound)) {
                                //std::cout << "setting upper bound in column " << r.k + shift_repair_col+1 << " to " << r.bound << std::endl;
                            }
                            else {
                                throw LPBConversionException("Upper bound could not be decreased.");
                            }
                        }
                        k = r.k + shift_repair_col + 1;
                        continue;
                    }
                }
                --k;
            }
            interval_t rootInterval =
                    information.getInterval(*(tree->getRoot()));
            lpb_coefficient_t degree =
                    chooseDegree(rootInterval.first, rootInterval.second);
            LPB *result = new LPB(degree, numVariables);
            VariableSetting *vs = tree->getVariableSetting();
            for (size_t i = 0; i < numVariables; i++) {
                size_t literalID = vs->operator[](i);
                result->addSummand(information.getCoefficient(i + 1),
                        Literal(literalID));
            }
            return result;
        }
    protected:
        /**
         * Construct an error message.
         *
         * The returned message says that it is not possible to
         * choose the degree in the interval (a, b].
         *
         * @param a First value of the interval.
         * @param b Second value of the interval.
         *
         * @return Error message.
         */
        inline std::string errorMessage(lpb_coefficient_t a,
                lpb_coefficient_t b) const {
            std::ostringstream error;
            error << "Can't choose degree in the interval (" <<
                    a << ", " << b << "].";
            return error.str();
        }

  //muc and mlc enforce a lower bound, nuc and nlc an upper bound. But it's only in the non-recursive
  //part where this matters. 
        RepairResult repair(RepairInput &muc, RepairInput &mlc,
                RepairInput &nuc, RepairInput &nlc,
                size_t column, size_t m) {
	  //          std::cout << "REPAIR " << std::endl;
            if (column == m)
                return RepairResult();
            interval_t muUI;
            interval_t muLI;
            interval_t mlUI;
            interval_t mlLI;
            interval_t nuUI;
            interval_t nuLI;
            interval_t nlUI;
            interval_t nlLI;
            std::pair<RepairInput, RepairInput> r1 = createResults(muc);
            //~ std::cout << "muc is final? " << muc.isFinal() << std::endl;
            //~ std::cout << "mlc is final? " << mlc.isFinal() << std::endl;
            //~ std::cout << "nuc is final? " << nuc.isFinal() << std::endl;
            //~ std::cout << "nlc is final? " << nlc.isFinal() << std::endl;
            //~ std::cout << "muU: isFinal: " << r1.first.isFinal() << ", trueType: " << r1.first.getContent().trueType << std::endl;
            //~ std::cout << "muL: isFinal: " << r1.second.isFinal() << ", trueType: " << r1.second.getContent().trueType << std::endl;
            
            std::pair<RepairInput, RepairInput> r2 = createResults(mlc);
            std::pair<RepairInput, RepairInput> r3 = createResults(nuc);
            std::pair<RepairInput, RepairInput> r4 = createResults(nlc);
            RepairInput muU = r1.first;
            RepairInput muL = r1.second;
            RepairInput mlU = r2.first;
            RepairInput mlL = r2.second;
            RepairInput nuU = r3.first;
            RepairInput nuL = r3.second;
            RepairInput nlU = r4.first;
            RepairInput nlL = r4.second;
            
            lpb_coefficient_t sum = information.getSumFrom(column + 2);
            if (muc.isFinal()
                    //~ BUGFIX second disjunct necessary due to cutTree option
		|| (muc.getContent().node->getUpperChild() == 0 || muc.getContent().node->getLowerChild() == 0))
	      {
                if (muc.getContent().trueType == 0) {
                    muUI = interval_t(sum, coeff_pos_infinity);
                    muLI = interval_t(sum, coeff_pos_infinity);
                } else {
                    muUI = interval_t(coeff_neg_infinity, 0);
                    muLI = interval_t(coeff_neg_infinity, 0);
                }
            } else {
                muUI = information.getInterval(*muc.getContent().node->getUpperChild());
                muLI = information.getInterval(*muc.getContent().node->getLowerChild());
            }
            //~ std::cout << "muUI: " << muUI << std::endl;
            //~ std::cout << "muLI: " << muLI << std::endl;

            if (mlc.isFinal()
                    //~ BUGFIX second disjunct necessary due to cutTree option
                    || (mlc.getContent().node->getUpperChild() == 0 || mlc.getContent().node->getLowerChild() == 0)) {
                if (mlc.getContent().trueType == 0) {
                    mlUI = interval_t(sum, coeff_pos_infinity);
                    mlLI = interval_t(sum, coeff_pos_infinity);
                } else {
                    mlUI = interval_t(coeff_neg_infinity, 0);
                    mlLI = interval_t(coeff_neg_infinity, 0);
                }
            } else {
                mlUI = information.getInterval(*mlc.getContent().node->getUpperChild());
                mlLI = information.getInterval(*mlc.getContent().node->getLowerChild());
            }
            //~ std::cout << "mlUI: " << mlUI << std::endl;
            //~ std::cout << "mlLI: " << mlLI << std::endl << std::endl;
            
            if (nuc.isFinal()
                    //~ BUGFIX second disjunct necessary due to cutTree option
                    || (nuc.getContent().node->getUpperChild() == 0 || nuc.getContent().node->getLowerChild() == 0)
                    ) {
                if (nuc.getContent().trueType == 0) {
                    nuUI = interval_t(sum, coeff_pos_infinity);
                    nuLI = interval_t(sum, coeff_pos_infinity);
                } else {
                    nuUI = interval_t(coeff_neg_infinity, 0);
                    nuLI = interval_t(coeff_neg_infinity, 0);
                }
            } else {
                nuUI = information.getInterval(*nuc.getContent().node->getUpperChild());
                nuLI = information.getInterval(*nuc.getContent().node->getLowerChild());
            }
            //std::cout << "nuUI: " << nuUI << std::endl;
            //~ std::cout << "nuLI: " << nuLI << std::endl;

            if (nlc.isFinal()
                    //~ BUGFIX second disjunct necessary due to cutTree option
                    || (nlc.getContent().node->getUpperChild() == 0 || nlc.getContent().node->getLowerChild() == 0)
                    ) {
                if (nlc.getContent().trueType == 0) {
                    nlUI = interval_t(sum, coeff_pos_infinity);
                    nlLI = interval_t(sum, coeff_pos_infinity);
                } else {
                    nlUI = interval_t(coeff_neg_infinity, 0);
                    nlLI = interval_t(coeff_neg_infinity, 0);
                }
            } else {
                nlUI = information.getInterval(*nlc.getContent().node->getUpperChild());
                nlLI = information.getInterval(*nlc.getContent().node->getLowerChild());
            }
            //std::cout << "nlUI: " << nlUI << std::endl;
            //~ std::cout << "nlLI: " << nlLI << std::endl << std::endl;

            bool need2Ms = false;
            bool need2Ns = false;

	    //Determine the bounds for sigma 
            RepairInput muPrimeSigma = muU;
            RepairInput nuPrimeSigma = muPrimeSigma;
            RepairInput mlPrimeSigma = mlU;
            RepairInput nlPrimeSigma = mlPrimeSigma;

	    //Determine the tightest lower bound for sigma
            lpb_coefficient_t l_sigma = muUI.first - mlUI.second;
            lpb_coefficient_t tmp;

            tmp = muLI.first - mlLI.second;
            if (tmp > l_sigma) {
                l_sigma = tmp;
                muPrimeSigma = muL;
                mlPrimeSigma = mlL;
            }

            tmp = nuUI.first - nlUI.second;
            if (tmp > l_sigma) {
                l_sigma = tmp;
                muPrimeSigma = nuU;
                mlPrimeSigma = nlU;
            }

            tmp = nuLI.first - nlLI.second;
            if (tmp > l_sigma) {
                l_sigma = tmp;
                muPrimeSigma = nuL;
                mlPrimeSigma = nlL;
            }

	    //Determine the tightest upper bound for sigma
            lpb_coefficient_t u_sigma = muUI.second - mlUI.first;

            tmp = muLI.second - mlLI.first;
            if (tmp < u_sigma) {
                u_sigma = tmp;
                nuPrimeSigma = muL;
                nlPrimeSigma = mlL;
            }

            tmp = nuUI.second - nlUI.first;
            if (tmp < u_sigma) {
                u_sigma = tmp;
                nuPrimeSigma = nuU;
                nlPrimeSigma = nlU;
            }

            tmp = nuLI.second - nlLI.first;
            if (tmp < u_sigma) {
                u_sigma = tmp;
                nuPrimeSigma = nuL;
                nlPrimeSigma = nlL;
            }

	    //std::cout << "u_sigma, l_sigma computed" << std::endl;
            //Conflict on sigma
            if (l_sigma >= u_sigma) {
	      //std::cout << "recursive call: case sigma" << std::endl;
	      //std::cout << muPrimeSigma.isFinal() << ", " << (bool)muPrimeSigma.getContent().trueType << std::endl;
                return repair(muPrimeSigma, mlPrimeSigma, nuPrimeSigma, nlPrimeSigma, column + 1, m);
            }

	    //Looking up the bounds for a_{k+2} is a way of deriving an additional constraint 
	    //for sigma + a_{k+2} and sigma - a_{k+2}.
            RepairInput muPrime_ak2(true);
            RepairInput nuPrime_ak2(true);
            RepairInput mlPrime_ak2(true);
            RepairInput nlPrime_ak2(true);

	    //Look up the bounds for sigma + a_{k+2}.
	    //Look up the tightest lower bound for a_{k+2}
	    lpb_coefficient_t l_ak2 = information.getComputedBound(column + 1).first;
	    //std::cout << "l_ak2 = " << l_ak2 << std::endl;


	    //Look up the tightest upper bound for a_{k+2}
	    lpb_coefficient_t u_ak2 = information.getComputedBound(column + 1).second;
	    //std::cout << "u_ak2 = " << u_ak2 << std::endl;

	    //Determine the bounds for sigma + a_{k+2}
            RepairInput muPrime = muU;
            RepairInput nuPrime = muPrime;
            RepairInput mlPrime = mlL;
            RepairInput nlPrime = mlPrime;

	    //Determine the tightest lower bound for sigma + a_{k+2}
            lpb_coefficient_t l_sigma_plus = muUI.first - mlLI.second;

            tmp = nuUI.first - nlLI.second;
            if (tmp > l_sigma_plus) {
                l_sigma_plus = tmp;
                muPrime = nuU;
                mlPrime = nlL;
            }

	    //JGS : additional bound by adding up the bounds for sigma and for a_{k+2}
	    tmp = l_sigma + l_ak2;
            if (tmp > l_sigma_plus) {
	        //std::cout << "l_sigma + l_ak2 gave tightest lower bound for sigma + a_{k+2}" << std::endl;
                l_sigma_plus = tmp;
                need2Ms = true;
            }

	    //Determine the tightest upper bound for sigma + a_{k+2}
            lpb_coefficient_t u_sigma_plus = muUI.second - mlLI.first;

            tmp = nuUI.second - nlLI.first;
            if (tmp < u_sigma_plus) {
                u_sigma_plus = tmp;
                nuPrime = nuU;
                nlPrime = nlL;
            }

	    //JGS : additional bound by adding up the bounds for sigma and for a_{k+2}
	    tmp = u_sigma + u_ak2;
            if (tmp < l_sigma_plus) {
	        //std::cout << "u_sigma + u_ak2 gave tightest upper bound for sigma + a_{k+2}" << std::endl;
                u_sigma_plus = tmp;
                need2Ns = true;
            }


	    //Conflict on sigma + a_{k+2}
	    //std::cout << "Checking for conflict on sigma + a_{k+2}" << std::endl;
            if (l_sigma_plus >= u_sigma_plus) {
	      //std::cout << "l_sigma_plus >= u_sigma_plus " << l_sigma_plus << " >= " << u_sigma_plus << std::endl;
	      if (need2Ms) {//need to find the precomputed nodes that impose the tightest constraint on a_{k+2}. 
		//BEGIN mic mac stuff
		SplittingTreeNode *maxCulprit = information.getCulprit(column + 1).second;//!!!

		//		RepairInput muPrime_ak2(true);
		//		RepairInput mlPrime_ak2(true);

		RepairInputContent c5;
		RepairInputContent c6;

		if (maxCulprit->getUpperChild() != 0) {
		  muPrime_ak2.setFinal(false);
		  c5.node = maxCulprit;
		} else {
		  if (maxCulprit->getDNF()->size() == 0) {
		    c5.trueType = 0; } else {
		    c5.trueType = 1;}
		}

		if (maxCulprit->getLowerParent()->getLowerChild()->getUpperChild() != 0) {
		  mlPrime_ak2.setFinal(false);
		  c6.node = maxCulprit->getLowerParent()->getLowerChild();
		} else {
		  if (maxCulprit->getLowerParent()->getLowerChild()->getDNF()->size() == 0) {
		    c6.trueType = 0; } else {
		    c6.trueType = 1;}
		}

		muPrime_ak2.setContent(c5);
		mlPrime_ak2.setContent(c6);
		//END mic mac stuff
	      }
	      if (need2Ns) {//need to find the precomputed nodes that impose the tightest constraint on a_{k+2}. 
		//BEGIN mic mac stuff
		SplittingTreeNode *minCulprit = information.getCulprit(column + 1).first;//!!!

		//		RepairInput nuPrime_ak2(true);
		//		RepairInput nlPrime_ak2(true);

		RepairInputContent c5;
		RepairInputContent c6;

		if (minCulprit->getUpperChild() != 0) {
		  nuPrime_ak2.setFinal(false);
		  c5.node = minCulprit;
		} else {
		  if (minCulprit->getDNF()->size() == 0) {
		    c5.trueType = 0; } else {
		    c5.trueType = 1;}
		}

		if (minCulprit->getLowerParent()->getLowerChild()->getUpperChild() != 0) {
		  nlPrime_ak2.setFinal(false);
		  c6.node = minCulprit->getLowerParent()->getLowerChild();
		} else {
		  if (minCulprit->getLowerParent()->getLowerChild()->getDNF()->size() == 0) {
		    c6.trueType = 0; } else {
		    c6.trueType = 1;}
		}

		nuPrime_ak2.setContent(c5);
		nlPrime_ak2.setContent(c6);
		//END mic mac stuff
	      }		
	      if (need2Ms) {
		if (need2Ns) {
		  //std::cout << "Need repair3, which is not implemented" << std::endl;//This is sometimes needed!!!. JGS: clean!
		  exit(0);
		  //		  return repair3(muPrimeSigma, mlPrimeSigma, muPrime_ak2, mlPrime_ak2, 
		  //		 nuPrimeSigma, nlPrimeSigma, nuPrime_ak2, nlPrime_ak2, column + 1, m);
	        } else {
		    //goofy: N gives the lower bound, M the upper bound
		    return repair2(nuPrime, nlPrime,
                                   muPrimeSigma, mlPrimeSigma, muPrime_ak2, mlPrime_ak2, column + 1, m, true, false);		    
	        } 
	      } else {
		if (need2Ns) {
		    return repair2(muPrime, mlPrime,
                                   nuPrimeSigma, nlPrimeSigma, nuPrime_ak2, nlPrime_ak2, column + 1, m, false, false);
	        } else {
		    //std::cout << "recursive call: case sigma + a_{k+2}" << std::endl;
                    return repair(muPrime, mlPrime, nuPrime, nlPrime, column + 1, m);
		}
	      }
	    }

	    //Determine the bounds for sigma - a_{k+2}
            need2Ms = false;
            need2Ns = false;
            muPrime = muL;
            nuPrime = muPrime;
            mlPrime = mlU;
            nlPrime = mlPrime;

	    //Determine the tightest lower bound for sigma - a_{k+2}
            lpb_coefficient_t l_sigma_minus = muLI.first - mlUI.second;

            tmp = nuLI.first - nlUI.second;
            if ( tmp > l_sigma_minus) {
                l_sigma_minus =  tmp;
                muPrime = nuL;
                mlPrime = nlU;
            }

	    //JGS : additional bound by taking the difference between the bounds for sigma and for a_{k+2}
	    tmp = l_sigma - u_ak2; //upper bound for a_{k+2} influences lower bound for sigma - a_{k+2}
            if (tmp > l_sigma_minus) {
	        //std::cout << "l_sigma - u_ak2 gave tightest lower bound for sigma - a_{k+2}" << std::endl;
                l_sigma_minus = tmp;
                need2Ms = true;
            }
	    //std::cout << "l_sigma_minus = " << l_sigma_minus << std::endl; 

	    //Determine the tightest upper bound for sigma - a_{k+2}
            lpb_coefficient_t u_sigma_minus = muLI.second - mlUI.first;

            tmp = nuLI.second - nlUI.first;
            if (tmp < u_sigma_minus) {
                u_sigma_minus = tmp;
                nuPrime = nuL;
                nlPrime = nlU;
            }

	    //JGS : additional bound by adding up the bounds for sigma and for a_{k+2}
	    tmp = u_sigma - l_ak2; //lower bound for a_{k+2} influences upper bound for sigma - a_{k+2}
	    //	    std::cout << "u_sigma = " << u_sigma << std::endl;
            if (tmp < u_sigma_minus) {
	        //std::cout << "u_sigma - l_ak2 gave tightest upper bound for sigma - a_{k+2}" << std::endl;
                u_sigma_minus = tmp;
                need2Ns = true;
            }
	    //std::cout << "u_sigma_minus = " << u_sigma_minus << std::endl;
 
	    //std::cout << "Checking for conflict on sigma - a_{k+2}" << std::endl; 
	    //Conflict on sigma - a_{k+2}
            if (l_sigma_minus >= u_sigma_minus) {
	      //std::cout << "Conflict on sigma - a_{k+2}" << std::endl;
	      if (need2Ns) {//need to find the precomputed nodes that impose the tightest constraint on a_{k+2}. 
		//BEGIN mic mac stuff: next 2 mic-mac sections are IDENTICAL to previous 2 mic-mac sections!
		SplittingTreeNode *maxCulprit = information.getCulprit(column + 1).second;//!!!
		//std::cout << "maxCulprit = " <<  information.getInterval(*maxCulprit) << std::endl;

		//		RepairInput muPrime_ak2(true);
		//		RepairInput mlPrime_ak2(true);

		RepairInputContent c5;
		RepairInputContent c6;

		if (maxCulprit->getUpperChild() != 0) {
		  muPrime_ak2.setFinal(false);
		  c5.node = maxCulprit;
		} else {
		  if (maxCulprit->getDNF()->size() == 0) {
		    c5.trueType = 0; } else {
		    c5.trueType = 1;}
		}

		if (maxCulprit->getLowerParent()->getLowerChild()->getUpperChild() != 0) {
		  mlPrime_ak2.setFinal(false);
		  c6.node = maxCulprit->getLowerParent()->getLowerChild();
		} else {
		  if (maxCulprit->getLowerParent()->getLowerChild()->getDNF()->size() == 0) {
		    c6.trueType = 0; } else {
		    //std::cout << "sibling of maxCulprit is true" << std::endl;
		    c6.trueType = 1;}
		}

		muPrime_ak2.setContent(c5);
		mlPrime_ak2.setContent(c6);
		//END mic mac stuff
	      }
	      if (need2Ms) {//need to find the precomputed nodes that impose the tightest constraint on a_{k+2}. 
		//BEGIN mic mac stuff
		SplittingTreeNode *minCulprit = information.getCulprit(column + 1).first;//!!!

		//		RepairInput nuPrime_ak2(true);
		//		RepairInput nlPrime_ak2(true);

		RepairInputContent c5;
		RepairInputContent c6;

		if (minCulprit->getUpperChild() != 0) {
		  nuPrime_ak2.setFinal(false);
		  c5.node = minCulprit;
		} else {
		  if (minCulprit->getDNF()->size() == 0) {
		    c5.trueType = 0; } else {
		    c5.trueType = 1;}
		}

		if (minCulprit->getLowerParent()->getLowerChild()->getUpperChild() != 0) {
		  nlPrime_ak2.setFinal(false);
		  c6.node = minCulprit->getLowerParent()->getLowerChild();
		} else {
		  if (minCulprit->getLowerParent()->getLowerChild()->getDNF()->size() == 0) {
		    c6.trueType = 0; } else {
		    c6.trueType = 1;}
		}

		nuPrime_ak2.setContent(c5);
		nlPrime_ak2.setContent(c6);
		//END mic mac stuff
	      }		
	      if (need2Ms) {
		if (need2Ns) {
		    return repair3(muPrimeSigma, mlPrimeSigma, muPrime_ak2, mlPrime_ak2, 
                                   nuPrimeSigma, nlPrimeSigma, nuPrime_ak2, nlPrime_ak2, column + 1, m);
	        } else {
		    //goofy: N gives the lower bound, M the upper bound
		    return repair2(nuPrime, nlPrime,
                                   muPrimeSigma, mlPrimeSigma, nuPrime_ak2, nlPrime_ak2, column + 1, m, true, true);		    
	        } 
	      } else {
		if (need2Ns) {
		    //std::cout << "Call to repair2 is imminent." << std::endl;
		    return repair2(muPrime, mlPrime,
                                   nuPrimeSigma, nlPrimeSigma, muPrime_ak2, mlPrime_ak2, column + 1, m, false, true);
	        } else {
		    //std::cout << "recursive call: case sigma - a_{k+2}" << std::endl;
                    return repair(muPrime, mlPrime, nuPrime, nlPrime, column + 1, m);
		}
	      }
	    }

            // non-recursive part
            lpb_coefficient_t ak2 = information.getCoefficient(column + 1);
	    //std::cout << "old coefficient: " << ak2 << std::endl;
            lpb_coefficient_t mus, mlb, nub, nls;
            lpb_coefficient_t leftA = 0;
            lpb_coefficient_t rightA = 0;
            if (muUI.first > (muLI.first + ak2)) {
                mus = muUI.first;
            }
            else {
                mus = muLI.first;
                ++leftA;
            }
	    //            std::cout << "mus: " << mus << std::endl;
	    //            std::cout << "leftA: " << leftA << std::endl;

            if (mlUI.second < (mlLI.second + ak2)) {
                mlb = mlUI.second;
            }
            else {
                mlb = mlLI.second;
                --leftA;
            }
            
	    //            std::cout << "mlb: " << mlb << std::endl;
	    //            std::cout << "leftA: " << leftA << std::endl;

            if (nuUI.second < (nuLI.second + ak2)) {
                nub = nuUI.second;
            }
            else {
                nub = nuLI.second;
                ++rightA;
            }
            
	    //            std::cout << "nub: " << nub << std::endl;
	    //            std::cout << "rightA: " << rightA << std::endl;

            if (nlUI.first > (nlLI.first + ak2)) {
                nls = nlUI.first;
            }
            else {
                nls = nlLI.first;
                --rightA;
            }
            
	    //            std::cout << "nls: " << nls << std::endl;
	    //            std::cout << "rightA: " << rightA << std::endl;

            const lpb_coefficient_t lhs = mus - mlb - nub + nls;
            const lpb_coefficient_t rhs = rightA - leftA;
            
	    //std::cout << "lhs: " << lhs << std::endl;
	    //std::cout << "rhs: " << rhs << std::endl;

            // TODO: correct?!
            assert(rhs != 0);
            const lpb_coefficient_t bound = lhs / rhs;

	    // if (rhs < 0) {//is upper bound
	    //   return RepairResult(true, !(rhs < 0), ceil(bound), column + 1); } else {
	    //   return RepairResult(true, !(rhs < 0), floor(bound), column + 1);
	    // }//observe that rounding the repair bounds means loosening them, which is ok 
	    //as long as we calculate with integer coefficients, although due to the doubling I have doubts
	    //about this... 

	    //std::cout << "Will return" << std::endl;
	    return RepairResult(true, !(rhs < 0), bound, column + 1);
	}


  //muc and mlc are a node pair at distance sigma as for repair.
  //(nuc,nlc), and (n2uc,n2lc) are 2 node pairs whose distances sum up (or, if "minus" is true,
  //differ by) to sigma.
  //Unlike for repair, muc and mlc do not necessarily enforce a lower bound. 
  //If muc and mlc enforce an upper bound and (nuc,nlc), and (n2uc,n2lc) a lower bound,
  //we call this "goofy"
        RepairResult repair2(RepairInput &muc, RepairInput &mlc,
                RepairInput &nuc, RepairInput &nlc,
                RepairInput &n2uc, RepairInput &n2lc,
		size_t column, size_t m, bool goofy, bool minus) {
	  //	    std::cout << "REPAIR2 " << goofy << minus << std::endl;
            if (column == m)
                return RepairResult();

	    //std::cout << "n2u: isFinal: " << n2uc.isFinal() << std::endl;
	    //std::cout << "n2l: isFinal: " << n2lc.isFinal() << std::endl;

            interval_t muUI;
            interval_t muLI;
            interval_t mlUI;
            interval_t mlLI;
            interval_t nuUI;
            interval_t nuLI;
            interval_t nlUI;
            interval_t nlLI;
            interval_t n2uUI;
            interval_t n2uLI;
            interval_t n2lUI;
            interval_t n2lLI;
            std::pair<RepairInput, RepairInput> r1 = createResults(muc);
            //~ std::cout << "muc is final? " << muc.isFinal() << std::endl;
            //~ std::cout << "muU: isFinal: " << r1.first.isFinal() << ", trueType: " << r1.first.getContent().trueType << std::endl;
            //~ std::cout << "muL: isFinal: " << r1.second.isFinal() << ", trueType: " << r1.second.getContent().trueType << std::endl;
            
            std::pair<RepairInput, RepairInput> r2 = createResults(mlc);
            std::pair<RepairInput, RepairInput> r3 = createResults(nuc);
            std::pair<RepairInput, RepairInput> r4 = createResults(nlc);
            std::pair<RepairInput, RepairInput> r5 = createResults(n2uc);
            std::pair<RepairInput, RepairInput> r6 = createResults(n2lc);
	    //std::cout << "n2uu: isFinal: " << r5.first.isFinal() << ", trueType: " << r5.first.getContent().trueType << std::endl;
            //std::cout << "n2ul: isFinal: " << r5.second.isFinal() << ", trueType: " << r5.second.getContent().trueType << std::endl;
	    //std::cout << "n2lu: isFinal: " << r6.first.isFinal() << ", trueType: " << r6.first.getContent().trueType << std::endl;
            //std::cout << "n2ll: isFinal: " << r6.second.isFinal() << ", trueType: " << r6.second.getContent().trueType << std::endl;

            RepairInput muU = r1.first;
            RepairInput muL = r1.second;
            RepairInput mlU = r2.first;
            RepairInput mlL = r2.second;
            RepairInput nuU = r3.first;
            RepairInput nuL = r3.second;
            RepairInput nlU = r4.first;
            RepairInput nlL = r4.second;
            RepairInput n2uU = r5.first;
            RepairInput n2uL = r5.second;
            RepairInput n2lU = r6.first;
            RepairInput n2lL = r6.second;
            
            lpb_coefficient_t sum = information.getSumFrom(column + 2);

            if (muc.isFinal()
                    //~ BUGFIX second disjunct necessary due to cutTree option
		|| (muc.getContent().node->getUpperChild() == 0 || muc.getContent().node->getLowerChild() == 0)) {
                if (muc.getContent().trueType == 0) {
                    muUI = interval_t(sum, coeff_pos_infinity);
                    muLI = interval_t(sum, coeff_pos_infinity);
                } else {
                    muUI = interval_t(coeff_neg_infinity, 0);
                    muLI = interval_t(coeff_neg_infinity, 0);
                }
            } else {
                muUI = information.getInterval(*muc.getContent().node->getUpperChild());
                muLI = information.getInterval(*muc.getContent().node->getLowerChild());
            }
	    //std::cout << "muUI: " << muUI << std::endl;
	    //std::cout << "muLI: " << muLI << std::endl;

            if (mlc.isFinal()
                    //~ BUGFIX second disjunct necessary due to cutTree option
                    || (mlc.getContent().node->getUpperChild() == 0 || mlc.getContent().node->getLowerChild() == 0)) {
                if (mlc.getContent().trueType == 0) {
                    mlUI = interval_t(sum, coeff_pos_infinity);
                    mlLI = interval_t(sum, coeff_pos_infinity);
                } else {
                    mlUI = interval_t(coeff_neg_infinity, 0);
                    mlLI = interval_t(coeff_neg_infinity, 0);
                }
            } else {
                mlUI = information.getInterval(*mlc.getContent().node->getUpperChild());
                mlLI = information.getInterval(*mlc.getContent().node->getLowerChild());
            }
	    //std::cout << "mlUI: " << mlUI << std::endl;
	    //std::cout << "mlLI: " << mlLI << std::endl << std::endl;

            if (nuc.isFinal()
                    //~ BUGFIX second disjunct necessary due to cutTree option
                    || (nuc.getContent().node->getUpperChild() == 0 || nuc.getContent().node->getLowerChild() == 0)) {
                if (nuc.getContent().trueType == 0) {
                    nuUI = interval_t(sum, coeff_pos_infinity);
                    nuLI = interval_t(sum, coeff_pos_infinity);
                } else {
                    nuUI = interval_t(coeff_neg_infinity, 0);
                    nuLI = interval_t(coeff_neg_infinity, 0);
                }
            } else {
                nuUI = information.getInterval(*nuc.getContent().node->getUpperChild());
                nuLI = information.getInterval(*nuc.getContent().node->getLowerChild());
            }
	    //std::cout << "nuUI: " << nuUI << std::endl;
	    //std::cout << "nuLI: " << nuLI << std::endl;

            if (nlc.isFinal()
                    //~ BUGFIX second disjunct necessary due to cutTree option
                    || (nlc.getContent().node->getUpperChild() == 0 || nlc.getContent().node->getLowerChild() == 0)) {
                if (nlc.getContent().trueType == 0) {
                    nlUI = interval_t(sum, coeff_pos_infinity);
                    nlLI = interval_t(sum, coeff_pos_infinity);
                } else {
                    nlUI = interval_t(coeff_neg_infinity, 0);
                    nlLI = interval_t(coeff_neg_infinity, 0);
                }
            } else {
                nlUI = information.getInterval(*nlc.getContent().node->getUpperChild());
                nlLI = information.getInterval(*nlc.getContent().node->getLowerChild());
            }
	    //std::cout << "nlUI: " << nlUI << std::endl;
	    //std::cout << "nlLI: " << nlLI << std::endl << std::endl;

            if (n2uc.isFinal()
                    //~ BUGFIX second disjunct necessary due to cutTree option
                    || (n2uc.getContent().node->getUpperChild() == 0 || n2uc.getContent().node->getLowerChild() == 0)) {
                if (n2uc.getContent().trueType == 0) {
                    n2uUI = interval_t(sum, coeff_pos_infinity);
                    n2uLI = interval_t(sum, coeff_pos_infinity);
                } else {
                    n2uUI = interval_t(coeff_neg_infinity, 0);
                    n2uLI = interval_t(coeff_neg_infinity, 0);
                }
            } else {
                n2uUI = information.getInterval(*n2uc.getContent().node->getUpperChild());
                n2uLI = information.getInterval(*n2uc.getContent().node->getLowerChild());
            }
	    //std::cout << "n2uUI: " << n2uUI << std::endl;
	    //std::cout << "n2uLI: " << n2uLI << std::endl;

            if (n2lc.isFinal()
                    //~ BUGFIX second disjunct necessary due to cutTree option
                    || (n2lc.getContent().node->getUpperChild() == 0 || n2lc.getContent().node->getLowerChild() == 0)) {
                if (n2lc.getContent().trueType == 0) {
                    n2lUI = interval_t(sum, coeff_pos_infinity);
                    n2lLI = interval_t(sum, coeff_pos_infinity);
                } else {
                    n2lUI = interval_t(coeff_neg_infinity, 0);
                    n2lLI = interval_t(coeff_neg_infinity, 0);
                }
            } else {
                n2lUI = information.getInterval(*n2lc.getContent().node->getUpperChild());
                n2lLI = information.getInterval(*n2lc.getContent().node->getLowerChild());
            }
	    //std::cout << "n2lUI: " << n2lUI << std::endl;
	    //std::cout << "n2lLI: " << n2lLI << std::endl << std::endl;

	    //Look up the bounds for sigma + a_{k+2}.
	    //Look up the tightest lower bound for a_{k+2}
	    lpb_coefficient_t l_ak2 = information.getComputedBound(column + 1).first;
	    //std::cout << "l_ak2 = " << l_ak2 << std::endl;


	    //Look up the tightest upper bound for a_{k+2}
	    lpb_coefficient_t u_ak2 = information.getComputedBound(column + 1).second;
	    //std::cout << "u_ak2 = " << u_ak2 << std::endl;

	    bool need2Ms = false;
	    bool need2Ns = false;

	    //Determine the bounds for a_rho2 coming from N2
	    //Determine the tightest lower bound for a_rho2
	    RepairInput muPrimeRho2 = n2uU; 
            RepairInput mlPrimeRho2 = n2lU; 
	    RepairInput nuPrimeRho2 = n2uU;
            RepairInput nlPrimeRho2 = n2lU;
	    lpb_coefficient_t l_rho2, u_rho2;
	    if (n2uUI.first - n2lUI.second > n2uLI.first - n2lLI.second) {
	      l_rho2 = n2uUI.first - n2lUI.second;
	      if (minus) {
		nuPrimeRho2 = n2uU;
		nlPrimeRho2 = n2lU; 
	      } else {
		muPrimeRho2 = n2uU;
		mlPrimeRho2 = n2lU; 
	      }
	    }
	    else {
	      l_rho2 = n2uLI.first - n2lLI.second;
	      if (minus) {
		nuPrimeRho2 = n2uL;
		nlPrimeRho2 = n2lL; 
	      } else {
		muPrimeRho2 = n2uL;
		mlPrimeRho2 = n2lL; 
	      }
	    }
	    //std::cout << "l_rho2 = " << l_rho2 << std::endl;
	      
	    //Determine the tightest upper bound for a_rho2
	    if (n2uUI.second - n2lUI.first < n2uLI.second - n2lLI.first) {
	      u_rho2 = n2uUI.second - n2lUI.first;
	      if (minus) {
		muPrimeRho2 = n2uU;
		mlPrimeRho2 = n2lU; 
	      } else {
		nuPrimeRho2 = n2uU;
		nlPrimeRho2 = n2lU; 
	      }
	    }
	    else {
	      u_rho2 = n2uLI.second - n2lLI.first;
	      if (minus) {
		muPrimeRho2 = n2uL;
		mlPrimeRho2 = n2lL;
	      } else {
		nuPrimeRho2 = n2uL;
		nlPrimeRho2 = n2lL;
	      }
	    }
	    //std::cout << "u_rho2 = " << u_rho2 << std::endl;

	    //Determine the bounds for a_rho2 - a_{k+2} coming from N2
	    lpb_coefficient_t l_rho2Minus =  n2uLI.first - n2lUI.second;
 	    lpb_coefficient_t u_rho2Minus =  n2uLI.second - n2lUI.first;
	    //std::cout << "l_rho2Minus = " << l_rho2Minus << std::endl;
	    //std::cout << "u_rho2Minus = " << u_rho2Minus << std::endl;

	    //Determine the bounds for a_rho2 + a_{k+2} coming from N2
	    lpb_coefficient_t l_rho2Plus =  n2uUI.first - n2lLI.second;
 	    lpb_coefficient_t u_rho2Plus =  n2uUI.second - n2lLI.first;
	    //std::cout << "l_rho2Plus = " << l_rho2Plus << std::endl;
	    //std::cout << "u_rho2Plus = " << u_rho2Plus << std::endl;

	    //Determine the bounds for sigma 
            RepairInput muPrimeSigma = muU;
            RepairInput nuPrimeSigma = muPrimeSigma;
            RepairInput mlPrimeSigma = mlU;
            RepairInput nlPrimeSigma = mlPrimeSigma;

	    //Determine the tightest lower bound for sigma
            lpb_coefficient_t l_sigma = muUI.first - mlUI.second;
	    //std::cout << "Initial l_sigma = " << l_sigma << std::endl;
            lpb_coefficient_t tmp;

            tmp = muLI.first - mlLI.second;
            if (tmp > l_sigma) {
                l_sigma = tmp;
		//std::cout << "l_sigma = " << l_sigma << std::endl;
                muPrimeSigma = muL;
                mlPrimeSigma = mlL;
            }
            
            if (minus) { tmp = (nuUI.first - nlUI.second) - u_rho2; }
	    else { tmp = (nuUI.first - nlUI.second) + l_rho2; }
            if (tmp > l_sigma) {
                l_sigma = tmp;
		//std::cout << "l_sigma = " << l_sigma << std::endl;
                muPrimeSigma = nuU;
                mlPrimeSigma = nlU;
		need2Ms = true;
            }

            if (minus) { tmp = (nuLI.first - nlLI.second) - u_rho2; }
	    else { tmp = (nuLI.first - nlLI.second) + l_rho2; }
            if (tmp > l_sigma) {
                l_sigma = tmp;
		//std::cout << "l_sigma = " << l_sigma << std::endl;
                muPrimeSigma = nuL;
                mlPrimeSigma = nlL;
		need2Ms = true;
            }

	    RepairInput muPrimeRho2Minus(true);
	    RepairInput mlPrimeRho2Minus(true);
	    bool needRho2Minus = false;
            if (minus) { //case sigma = (rho - a_{k+2}) - (rho2 - a_{k+2})
	      tmp = (nuLI.first - nlUI.second) - u_rho2Minus; 
	      if (tmp > l_sigma) {
                l_sigma = tmp;
		//std::cout << "case sigma = (rho - a_{k+2}) - (rho2 - a_{k+2}). l_sigma = " << l_sigma << std::endl;
                muPrimeSigma = nuL;
                mlPrimeSigma = nlU;
		muPrimeRho2Minus = n2uL;
		mlPrimeRho2Minus = n2lU;
		needRho2Minus = true;
		need2Ms = true;
	      }
	    } else { //case sigma = (rho + a_{k+2}) + (rho2 - a_{k+2})
	      tmp = (nuUI.first - nlLI.second) + l_rho2Minus;
	      if (tmp > l_sigma) {
                l_sigma = tmp;
		//std::cout << "case sigma = (rho + a_{k+2}) + (rho2 - a_{k+2}) l_sigma = " << l_sigma << std::endl;
                muPrimeSigma = nuU;
                mlPrimeSigma = nlL;
		muPrimeRho2Minus = n2uL;
		mlPrimeRho2Minus = n2lU;
		needRho2Minus = true;
		need2Ms = true;
	      }
            }

	    RepairInput muPrimeRho2Plus(true);
	    RepairInput mlPrimeRho2Plus(true);
	    bool needRho2Plus = false;
            if (minus) { //case sigma = (rho + a_{k+2}) - (rho2 + a_{k+2})
	      tmp = (nuUI.first - nlLI.second) - u_rho2Plus; 
	      if (tmp > l_sigma) {
                l_sigma = tmp;
		//std::cout << "case sigma = (rho + a_{k+2}) - (rho2 + a_{k+2}). l_sigma = " << l_sigma << std::endl;
                muPrimeSigma = nuU;
                mlPrimeSigma = nlL;
		muPrimeRho2Plus = n2uU;
		mlPrimeRho2Plus = n2lL;
		needRho2Plus = true;
		need2Ms = true;
	      }
	    } else { //case sigma = (rho - a_{k+2}) + (rho2 + a_{k+2})
	      tmp = (nuLI.first - nlUI.second) + l_rho2Plus;
	      if (tmp > l_sigma) {
                l_sigma = tmp;
		//std::cout << "case sigma = (rho - a_{k+2}) + (rho2 + a_{k+2}). l_sigma = " << l_sigma << std::endl;
                muPrimeSigma = nuL;
                mlPrimeSigma = nlU;
		muPrimeRho2Plus = n2uU;
		mlPrimeRho2Plus = n2lL;
		needRho2Plus = true;
		need2Ms = true;
	      }
            }
	    //std::cout << "Final l_sigma = " << l_sigma << std::endl;

	    //Determine the tightest upper bound for sigma
            lpb_coefficient_t u_sigma = muUI.second - mlUI.first;
	    //std::cout << "Initial u_sigma = " << u_sigma << std::endl;

            tmp = muLI.second - mlLI.first;
            if (tmp < u_sigma) {
                u_sigma = tmp;
                nuPrimeSigma = muL;
                nlPrimeSigma = mlL;
            }

            if (minus) { tmp = (nuUI.second - nlUI.first) - l_rho2; }
	    else { tmp = (nuUI.second - nlUI.first) + u_rho2; }
            if (tmp < u_sigma) {
                u_sigma = tmp;
                nuPrimeSigma = nuU;
                nlPrimeSigma = nlU;
		need2Ns = true;
            }

            if (minus) { tmp = (nuLI.second - nlLI.first) - l_rho2; }
	    else { tmp = (nuLI.second - nlLI.first) + u_rho2; }
            if (tmp < u_sigma) {
                u_sigma = tmp;
                nuPrimeSigma = nuL;
                nlPrimeSigma = nlL;
		need2Ns = true;
            }

	    RepairInput nuPrimeRho2Minus(true);
	    RepairInput nlPrimeRho2Minus(true);
            if (minus) { //case sigma = (rho - a_{k+2}) - (rho2 - a_{k+2})
	      tmp = (nuLI.second - nlUI.first) - l_rho2Minus; 
	      if (tmp < u_sigma) {
                u_sigma = tmp;
		//std::cout << "case sigma = (rho - a_{k+2}) - (rho2 - a_{k+2}). u_sigma = " << u_sigma << std::endl;
                nuPrimeSigma = nuL;
                nlPrimeSigma = nlU;
		nuPrimeRho2Minus = n2uL;
		nlPrimeRho2Minus = n2lU;
		needRho2Minus = true;
		need2Ns = true;
	      }
	    } else { //case sigma = (rho + a_{k+2}) + (rho2 - a_{k+2})
	      tmp = (nuUI.second - nlLI.first) + u_rho2Minus;
	      if (tmp < u_sigma) {
                u_sigma = tmp;
		//std::cout << "case sigma = (rho + a_{k+2}) + (rho2 - a_{k+2}). u_sigma = " << u_sigma << std::endl;
                nuPrimeSigma = nuU;
                nlPrimeSigma = nlL;
		nuPrimeRho2Minus = n2uL;
		nlPrimeRho2Minus = n2lU;
		needRho2Minus = true;
		need2Ns = true;
	      }
            }

	    RepairInput nuPrimeRho2Plus(true);
	    RepairInput nlPrimeRho2Plus(true);
            if (minus) { //case sigma = (rho + a_{k+2}) - (rho2 + a_{k+2})
	      tmp = (nuUI.second - nlLI.first) - l_rho2Plus; 
	      if (tmp < u_sigma) {
                u_sigma = tmp;
		//std::cout << "case sigma = (rho + a_{k+2}) - (rho2 + a_{k+2}). u_sigma = " << u_sigma << std::endl;
                nuPrimeSigma = nuU;
                nlPrimeSigma = nlL;
		nuPrimeRho2Plus = n2uU;
		nlPrimeRho2Plus = n2lL;
		needRho2Plus = true;
		need2Ns = true;
	      }
	    } else { //case sigma = (rho - a_{k+2}) + (rho2 + a_{k+2})
	      tmp = (nuLI.second - nlUI.first) + u_rho2Plus;
	      if (tmp < u_sigma) {
                u_sigma = tmp;
		//std::cout << "u_sigma = " << u_sigma << std::endl;
                nuPrimeSigma = nuL;
                nlPrimeSigma = nlU;
		nuPrimeRho2Plus = n2uU;
		nlPrimeRho2Plus = n2lL;
		needRho2Plus = true;
		need2Ns = true;
	      }
            }
	    //std::cout << "Final u_sigma = " << u_sigma << std::endl;

	    //Conflict on sigma
            if (l_sigma >= u_sigma) {
	      if (needRho2Plus) {
		muPrimeRho2 = muPrimeRho2Plus;
		mlPrimeRho2 = mlPrimeRho2Plus;
		nuPrimeRho2 = nuPrimeRho2Plus;
		nlPrimeRho2 = nlPrimeRho2Plus;
	      } else {
		if (needRho2Minus) {
		  muPrimeRho2 = muPrimeRho2Minus;
		  mlPrimeRho2 = mlPrimeRho2Minus;
		  nuPrimeRho2 = nuPrimeRho2Minus;
		  nlPrimeRho2 = nlPrimeRho2Minus;
		}
	      }
	      if (need2Ms) {
		if (need2Ns) {
		    return repair3(muPrimeSigma, mlPrimeSigma, muPrimeRho2, mlPrimeRho2, 
                                   nuPrimeSigma, nlPrimeSigma, nuPrimeRho2, nlPrimeRho2, column + 1, m);
	        } else {
		    //goofy: N gives the lower bound, M the upper bound
		    return repair2(nuPrimeSigma, nlPrimeSigma,
                                   muPrimeSigma, mlPrimeSigma, muPrimeRho2, mlPrimeRho2, column + 1, m, true, minus);		    
	        } 
	      } else {
		if (need2Ns) {
		    return repair2(muPrimeSigma, mlPrimeSigma,
                                   nuPrimeSigma, nlPrimeSigma, nuPrimeRho2, nlPrimeRho2, column + 1, m, false, minus);
	        } else {
		    //std::cout << "recursive call: case 2" << std::endl;
                    return repair(muPrimeSigma, mlPrimeSigma, nuPrimeSigma, nlPrimeSigma, column + 1, m);
		}
	      }
            } else {
	      //std::cout << "No conflict on sigma" << std::endl; 
	    }

	    //Determine the bounds for sigma + a_{k+2}
            RepairInput muPrime = muU;
            RepairInput nuPrime = muPrime;
            RepairInput mlPrime = mlL;
            RepairInput nlPrime = mlPrime;

	    //NEW
	    needRho2Plus = false;
	    needRho2Minus = false;

	    //Determine the tightest lower bound for sigma + a_{k+2}
            lpb_coefficient_t l_sigma_plus = muUI.first - mlLI.second;
	    if (minus) { tmp = (nuUI.first - nlLI.second) - u_rho2; }
	    else { tmp = (nuUI.first - nlLI.second) + l_rho2; }
            if (tmp > l_sigma_plus) {
                l_sigma_plus = tmp;
                muPrime = nuU;
                mlPrime = nlL;
		need2Ms = true;
            }
	    //NEW
            if (minus) { //case sigma + a_{k+2} = rho - (rho2 - a_{k+2})
	      tmp = (nuUI.first - nlUI.second) - u_rho2Minus; 
	      if (tmp > l_sigma_plus) {
                l_sigma_plus = tmp;
		//std::cout << "case sigma + a_{k+2} = rho - (rho2 - a_{k+2}). l_sigma_plus = " << l_sigma_plus << std::endl;
                muPrime = nuU;
                mlPrime = nlU;
		muPrimeRho2Minus = n2uL;
		mlPrimeRho2Minus = n2lU;
		needRho2Minus = true;
		need2Ms = true;
	      }
	      tmp = (nuLI.first - nlLI.second) - u_rho2Minus; 
	      if (tmp > l_sigma_plus) {
                l_sigma_plus = tmp;
		//std::cout << "case sigma + a_{k+2} = rho - (rho2 - a_{k+2}). l_sigma_plus = " << l_sigma_plus << std::endl;
                muPrime = nuL;
                mlPrime = nlL;
		muPrimeRho2Minus = n2uL;
		mlPrimeRho2Minus = n2lU;
		needRho2Minus = true;
		need2Ms = true;
	      }
	    } else { //case sigma + a_{k+2} = rho + (rho2 + a_{k+2})
	      tmp = (nuUI.first - nlUI.second) + l_rho2Plus;
	      if (tmp > l_sigma_plus) {
                l_sigma_plus = tmp;
		//std::cout << "case sigma + a_{k+2} = rho + (rho2 + a_{k+2}) l_sigma_plus = " << l_sigma_plus << std::endl;
                muPrime = nuU;
                mlPrime = nlU;
		muPrimeRho2Plus = n2uU;
		mlPrimeRho2Plus = n2lL;
		needRho2Plus = true;
		need2Ms = true;
	      }
	      tmp = (nuLI.first - nlLI.second) + l_rho2Plus;
	      if (tmp > l_sigma_plus) {
                l_sigma_plus = tmp;
		//std::cout << "case sigma + a_{k+2} = rho + (rho2 + a_{k+2}) l_sigma_plus = " << l_sigma_plus << std::endl;
                muPrime = nuL;
                mlPrime = nlL;
		muPrimeRho2Plus = n2uU;
		mlPrimeRho2Plus = n2lL;
		needRho2Plus = true;
		need2Ms = true;
	      }
            }
	    //END NEW
	    //std::cout << "l_sigma_plus = " << l_sigma_plus << std::endl;


	    //Determine the tightest upper bound for sigma + a_{k+2}
            lpb_coefficient_t u_sigma_plus = muUI.second - mlLI.first;
	    if (minus) { tmp = (nuUI.second - nlLI.first) - l_rho2; }
	    else { tmp = (nuUI.second - nlLI.first) + u_rho2; }
            if (tmp < u_sigma_plus) {
                u_sigma_plus = tmp;
                nuPrime = nuU;
                nlPrime = nlL;
		need2Ns = true;
            }
	    //NEW
            if (minus) { //case sigma + a_{k+2} = rho - (rho2 - a_{k+2})
	      tmp = (nuUI.second - nlUI.first) - l_rho2Minus; 
	      if (tmp < u_sigma_plus) {
                u_sigma_plus = tmp;
		//std::cout << "case sigma + a_{k+2} = rho - (rho2 - a_{k+2}). u_sigma_plus = " << u_sigma_plus << std::endl;
                nuPrime = nuU;
                nlPrime = nlU;
		nuPrimeRho2Minus = n2uL;
		nlPrimeRho2Minus = n2lU;
		needRho2Minus = true;
		need2Ns = true;
	      }
	      tmp = (nuLI.second - nlLI.first) - l_rho2Minus; 
	      if (tmp < u_sigma_plus) {
                u_sigma_plus = tmp;
		//std::cout << "case sigma + a_{k+2} = rho - (rho2 - a_{k+2}). u_sigma_plus = " << u_sigma_plus << std::endl;
                nuPrime = nuL;
                nlPrime = nlL;
		nuPrimeRho2Minus = n2uL;
		nlPrimeRho2Minus = n2lU;
		needRho2Minus = true;
		need2Ns = true;
	      }
	    } else { //case sigma + a_{k+2} = rho + (rho2 + a_{k+2})
	      tmp = (nuUI.second - nlUI.first) + u_rho2Plus;
	      if (tmp < u_sigma_plus) {
                u_sigma_plus = tmp;
		//std::cout << "case sigma + a_{k+2} = rho + (rho2 + a_{k+2}). u_sigma_plus = " << u_sigma_plus << std::endl;
                nuPrime = nuU;
                nlPrime = nlU;
		nuPrimeRho2Plus = n2uU;
		nlPrimeRho2Plus = n2lL;
		needRho2Plus = true;
		need2Ns = true;
	      }
	      tmp = (nuLI.second - nlLI.first) + u_rho2Plus;
	      if (tmp < u_sigma_plus) {
                u_sigma_plus = tmp;
		//std::cout << "case sigma + a_{k+2} = rho + (rho2 + a_{k+2}). u_sigma_plus = " << u_sigma_plus << std::endl;
                nuPrime = nuL;
                nlPrime = nlL;
		nuPrimeRho2Plus = n2uU;
		nlPrimeRho2Plus = n2lL;
		needRho2Plus = true;
		need2Ns = true;
	      }
            }
	    //END NEW
	    //std::cout << "u_sigma_plus = " << u_sigma_plus << std::endl;

	    //Conflict on sigma + a_{k+2}
            if (l_sigma_plus >= u_sigma_plus) {
	      if (needRho2Plus) {
		//std::cout << "assign PrimeRho2 = PrimeRho2Plus" << std::endl;
		muPrimeRho2 = muPrimeRho2Plus;
		mlPrimeRho2 = mlPrimeRho2Plus;
		nuPrimeRho2 = nuPrimeRho2Plus;
		nlPrimeRho2 = nlPrimeRho2Plus;
	      } else {
		if (needRho2Minus) {
		  //std::cout << "assign PrimeRho2 = PrimeRho2Minus" << std::endl;
		  muPrimeRho2 = muPrimeRho2Minus;
		  mlPrimeRho2 = mlPrimeRho2Minus;
		  nuPrimeRho2 = nuPrimeRho2Minus;
		  nlPrimeRho2 = nlPrimeRho2Minus;
		}
	      }
	      if (need2Ms) {
		if (need2Ns) {
		    return repair3(muPrime, mlPrime, muPrimeRho2, mlPrimeRho2, 
                                   nuPrime, nlPrime, nuPrimeRho2, nlPrimeRho2, column + 1, m);
	          } else {
		    //goofy: N gives the lower bound, M the upper bound
		    return repair2(nuPrime, nlPrime,
                                   muPrime, mlPrime, muPrimeRho2, mlPrimeRho2, column + 1, m, true, minus);		    
		  } 
	      } else {
		if (need2Ns) {
		    return repair2(muPrime, mlPrime,
                                   nuPrime, nlPrime, nuPrimeRho2, nlPrimeRho2, column + 1, m, false, minus);
	          } else {
		    //std::cout << "recursive call in repair2: case 2" << std::endl;
                    return repair(muPrime, mlPrime, nuPrime, nlPrime, column + 1, m);
		}
	      }
	    } else {
	      //std::cout << "No conflict on sigma + a_{k+2}" << std::endl; 
	    }

	    //Determine the bounds for sigma - a_{k+2}
            muPrime = muL;
            nuPrime = muPrime;
            mlPrime = mlU;
            nlPrime = mlPrime;

	    //NEW
	    needRho2Plus = false;
	    needRho2Minus = false;

	    //Determine the tightest lower bound for sigma - a_{k+2}
            lpb_coefficient_t l_sigma_minus = muLI.first - mlUI.second;
	    if (minus) { tmp = (nuLI.first - nlUI.second) - u_rho2; }
            else { tmp = (nuLI.first - nlUI.second) + l_rho2; }
            if ( tmp > l_sigma_minus) {
	      l_sigma_minus =  tmp;
	      //std::cout << "l_sigma_minus updated = " << l_sigma_minus << std::endl;
	      muPrime = nuL;
	      mlPrime = nlU;
	      need2Ms = true;
            }

	    //NEW
            if (minus) { //case sigma - a_{k+2} = rho - (rho2 + a_{k+2})
	      tmp = (nuUI.first - nlUI.second) - u_rho2Plus; 
	      if (tmp > l_sigma_minus) {
                l_sigma_minus = tmp;
		//std::cout << "case sigma - a_{k+2} = rho - (rho2 + a_{k+2}). l_sigma_minus = " << l_sigma_minus << std::endl;
                muPrime = nuU;
                mlPrime = nlU;
		muPrimeRho2Plus = n2uU;
		mlPrimeRho2Plus = n2lL;
		needRho2Plus = true;
		need2Ms = true;
	      }
	      tmp = (nuLI.first - nlLI.second) - u_rho2Plus; 
	      if (tmp > l_sigma_minus) {
                l_sigma_minus = tmp;
		//std::cout << "case sigma - a_{k+2} = rho - (rho2 + a_{k+2}). l_sigma_minus = " << l_sigma_minus << std::endl;
                muPrime = nuL;
                mlPrime = nlL;
		muPrimeRho2Plus = n2uU;
		mlPrimeRho2Plus = n2lL;
		needRho2Plus = true;
		need2Ms = true;
	      }
	    } else { //case sigma - a_{k+2} = rho + (rho2 - a_{k+2})
	      tmp = (nuUI.first - nlUI.second) + l_rho2Minus;
	      //std::cout << "(nuUI.first - nlUI.second) = " << (nuUI.first - nlUI.second) <<
	      //". l_rho2Minus = " << l_rho2Minus << std::endl;
	      if (tmp > l_sigma_minus) {
                l_sigma_minus = tmp;
		//std::cout << "case sigma - a_{k+2} = rho + (rho2 - a_{k+2}) l_sigma_minus = " << l_sigma_minus << std::endl;
                muPrime = nuU;
                mlPrime = nlU;
		muPrimeRho2Minus = n2uL;
		mlPrimeRho2Minus = n2lU;
		needRho2Minus = true;
		need2Ms = true;
	      }
	      tmp = (nuLI.first - nlLI.second) + l_rho2Minus;
	      if (tmp > l_sigma_minus) {
                l_sigma_minus = tmp;
		//std::cout << "case sigma - a_{k+2} = rho + (rho2 - a_{k+2}) l_sigma_minus = " << l_sigma_minus << std::endl;
                muPrime = nuL;
                mlPrime = nlL;
		muPrimeRho2Minus = n2uL;
		mlPrimeRho2Minus = n2lU;
		needRho2Minus = true;
		need2Ms = true;
	      }
            }
	    //END NEW
	    //std::cout << "l_sigma_minus = " << l_sigma_minus << std::endl;

	    //Determine the tightest upper bound for sigma - a_{k+2}
            lpb_coefficient_t u_sigma_minus = muLI.second - mlUI.first;
	    if (minus) { tmp = (nuLI.second - nlUI.first) - l_rho2; }
	    else { tmp = (nuLI.second - nlUI.first) + u_rho2; }
            if (tmp < u_sigma_minus) {
	      u_sigma_minus = tmp;
	      //std::cout << "u_sigma_minus updated = " << u_sigma_minus << std::endl;
	      nuPrime = nuL;
	      nlPrime = nlU;
	      need2Ns = true;
            }

	    //NEW
            if (minus) { //case sigma - a_{k+2} = rho - (rho2 + a_{k+2})
	      tmp = (nuUI.second - nlUI.first) - l_rho2Plus; 
	      if (tmp < u_sigma_minus) {
                u_sigma_minus = tmp;
		//std::cout << "case sigma - a_{k+2} = rho - (rho2 + a_{k+2}). u_sigma_minus = " << u_sigma_minus << std::endl;
                nuPrime = nuU;
                nlPrime = nlU;
		nuPrimeRho2Plus = n2uU;
		nlPrimeRho2Plus = n2lL;
		needRho2Plus = true;
		need2Ns = true;
	      }
	      tmp = (nuLI.second - nlLI.first) - l_rho2Plus; 
	      if (tmp < u_sigma_minus) {
                u_sigma_minus = tmp;
		//std::cout << "case sigma - a_{k+2} = rho - (rho2 + a_{k+2}). u_sigma_minus = " << u_sigma_minus << std::endl;
                nuPrime = nuL;
                nlPrime = nlL;
		nuPrimeRho2Plus = n2uU;
		nlPrimeRho2Plus = n2lL;
		needRho2Plus = true;
		need2Ns = true;
	      }
	    } else { //case sigma - a_{k+2} = rho + (rho2 - a_{k+2})
	      tmp = (nuUI.second - nlUI.first) + u_rho2Minus;
	      if (tmp < u_sigma_minus) {
                u_sigma_minus = tmp;
		//std::cout << "case sigma - a_{k+2} = rho + (rho2 - a_{k+2}). u_sigma_minus = " << u_sigma_minus << std::endl;
                nuPrime = nuU;
                nlPrime = nlU;
		nuPrimeRho2Minus = n2uL;
		nlPrimeRho2Minus = n2lU;
		needRho2Minus = true;
		need2Ns = true;
	      }
	      tmp = (nuLI.second - nlLI.first) + u_rho2Minus;
	      if (tmp < u_sigma_minus) {
                u_sigma_minus = tmp;
		//std::cout << "case sigma - a_{k+2} = rho + (rho2 - a_{k+2}). u_sigma_minus = " << u_sigma_minus << std::endl;
                nuPrime = nuL;
                nlPrime = nlL;
		nuPrimeRho2Minus = n2uL;
		nlPrimeRho2Minus = n2lU;
		needRho2Minus = true;
		need2Ns = true;
	      }
            }
	    //END NEW
	    //std::cout << "u_sigma_minus = " << u_sigma_minus << std::endl;

	    //Conflict on sigma - a_{k+2}
            if (l_sigma_minus >= u_sigma_minus) {
	      if (needRho2Plus) {
		//std::cout << "assign PrimeRho2 = PrimeRho2Plus" << std::endl;
		muPrimeRho2 = muPrimeRho2Plus;
		mlPrimeRho2 = mlPrimeRho2Plus;
		nuPrimeRho2 = nuPrimeRho2Plus;
		nlPrimeRho2 = nlPrimeRho2Plus;
	      } else {
		if (needRho2Minus) {
		  //std::cout << "assign PrimeRho2 = PrimeRho2Minus" << std::endl;
		  muPrimeRho2 = muPrimeRho2Minus;
		  mlPrimeRho2 = mlPrimeRho2Minus;
		  nuPrimeRho2 = nuPrimeRho2Minus;
		  nlPrimeRho2 = nlPrimeRho2Minus;
		}
	      }
	      if (need2Ms) {
		if (need2Ns) {
		    return repair3(muPrime, mlPrime, muPrimeRho2, mlPrimeRho2, 
                                   nuPrime, nlPrime, nuPrimeRho2, nlPrimeRho2, column + 1, m);
	          } else {
		    //goofy: N gives the lower bound, M the upper bound
		    return repair2(nuPrime, nlPrime,
                                   muPrime, mlPrime, muPrimeRho2, mlPrimeRho2, column + 1, m, true, minus);		    
		  } 
	      } else {
		if (need2Ns) {
		    return repair2(muPrime, mlPrime,
                                   nuPrime, nlPrime, nuPrimeRho2, nlPrimeRho2, column + 1, m, false, minus);
	          } else {
		    //std::cout << "recursive call in repair2: case 2" << std::endl;
                    return repair(muPrime, mlPrime, nuPrime, nlPrime, column + 1, m);
		}
	      }
            }
	    //std::cout << "No conflict on sigma - a_{k+2}" << std::endl;


            // non-recursive part
            lpb_coefficient_t ak2 = information.getCoefficient(column + 1);
	    //std::cout << "old coefficient: " << ak2 << std::endl;
	    //M lower bound, N upper bound
            lpb_coefficient_t mus, mlb, nub, nls, n2diff;
            lpb_coefficient_t leftA = 0;
            lpb_coefficient_t rightA = 0;
            lpb_coefficient_t n2A = 0;


	    if ( goofy == minus ) { //the combinatorics make that we need this condition
	      if (n2uUI.second < (n2uLI.second + ak2)) {
		n2diff = n2uUI.second;
	      } else {
		n2diff = n2uLI.second;
		++n2A;
	      }
	      //std::cout << "n2diff: " << n2diff << std::endl;
	      //std::cout << "n2A: " << n2A << std::endl;

	      if (n2lUI.first > (n2lLI.first + ak2)) {
		n2diff -= n2lUI.first;
	      } else {
		n2diff -= n2lLI.first;
		--n2A;
	      }
	      //std::cout << "n2diff: " << n2diff << std::endl;
	      //std::cout << "n2A: " << n2A << std::endl;

	      if (goofy) {
		n2diff *= -1;
		n2A *= -1;
	      }
	    } else {
	      if (n2lUI.second < (n2lLI.second + ak2)) {
		n2diff = n2lUI.second;
	      } else {
		n2diff = n2lLI.second;
		++n2A;
	      }
	      //std::cout << "n2diff: " << n2diff << std::endl;
	      //std::cout << "n2A: " << n2A << std::endl;

	      if (n2uUI.first > (n2uLI.first + ak2)) {
		  n2diff -= n2uUI.first;
	      } else {
		n2diff -= n2uLI.first;
		--n2A;
	      }
	      //std::cout << "n2diff: " << n2diff << std::endl;
	      //std::cout << "n2A: " << n2A << std::endl;

	      if (goofy) {
		n2diff *= -1;
		n2A *= -1;
	      }
	    }

	    if (goofy) {
	      mus = n2diff;
	      leftA = n2A;

	      if (nuUI.first > (nuLI.first + ak2)) {
		mus += nuUI.first;
	      } else {
		mus += nuLI.first;
		++leftA;
	      }
	      //std::cout << "mus: " << mus << std::endl;
	      //std::cout << "leftA: " << leftA << std::endl;

	      if (nlUI.second < (nlLI.second + ak2)) {
		mlb = nlUI.second;
	      } else {
		mlb = nlLI.second;
		--leftA;
	      }
	      //std::cout << "mlb: " << mlb << std::endl;
	      //std::cout << "leftA: " << leftA << std::endl;

	      if (muUI.second < (muLI.second + ak2)) {
		nub = muUI.second;
	      } else {
		nub = muLI.second;
		++rightA;
	      }
            
	      //std::cout << "nub: " << nub << std::endl;
	      //std::cout << "rightA: " << rightA << std::endl;

	      if (mlUI.first > (mlLI.first + ak2)) {
		nls = mlUI.first;
	      } else {
		nls = mlLI.first;
		--rightA;
	      }
            
	      //std::cout << "nls: " << nls << std::endl;
	      //std::cout << "rightA: " << rightA << std::endl;
	    } else {
	      nub = n2diff;
	      rightA = n2A;

		if (muUI.first > (muLI.first + ak2)) {
		  mus = muUI.first;
		} else {
		  mus = muLI.first;
		  ++leftA;
		}
		//std::cout << "mus: " << mus << std::endl;
		//std::cout << "leftA: " << leftA << std::endl;

		if (mlUI.second < (mlLI.second + ak2)) {
		  mlb = mlUI.second;
		} else {
		  mlb = mlLI.second;
		  --leftA;
		}
            
		//std::cout << "mlb: " << mlb << std::endl;
		//std::cout << "leftA: " << leftA << std::endl;

		if (nuUI.second < (nuLI.second + ak2)) {
		  nub += nuUI.second;
		} else {
		  nub += nuLI.second;
		  ++rightA;
		}
            
		//std::cout << "nub: " << nub << std::endl;
		//std::cout << "rightA: " << rightA << std::endl;

		if (nlUI.first > (nlLI.first + ak2)) {
		  nls = nlUI.first;
		} else {
		  nls = nlLI.first;
		  --rightA;
		}
            
		//std::cout << "nls: " << nls << std::endl;
		//std::cout << "rightA: " << rightA << std::endl;
	    }


	    const lpb_coefficient_t lhs = mus - mlb - nub + nls;
	    const lpb_coefficient_t rhs = rightA - leftA;
            
	    //std::cout << "lhs: " << lhs << std::endl;
	    //std::cout << "rhs: " << rhs << std::endl;

            // TODO: correct?!
            assert(rhs != 0);
            const lpb_coefficient_t bound = lhs / rhs;

	    // if (rhs < 0) {//is upper bound
	    //   return RepairResult(true, !(rhs < 0), ceil(bound), column + 1); } else {
	    //   return RepairResult(true, !(rhs < 0), floor(bound), column + 1);
	    // }

	    return RepairResult(true, !(rhs < 0), bound, column + 1);
	}

  
        RepairResult repair3(RepairInput &muc, RepairInput &mlc, RepairInput &m2uc, RepairInput &m2lc,
                RepairInput &nuc, RepairInput &nlc, RepairInput &n2uc, RepairInput &n2lc,
                size_t column, size_t m) {
	  std::cout << "repair3 not implemented yet." << std::endl;
	  return RepairResult();
	}

  //muc and mlc are a node pair at distance sigma as for repair.
  //Unlike for repair, muc and mlc do not necessarily enforce a lower bound. 
  //If muc and mlc enforce an upper bound, we call this "goofy"
        RepairResult monorepair(RepairInput &muc, RepairInput &mlc,
				size_t column, size_t m, lpb_coefficient_t oldbound, bool goofy) {
	  //	    std::cout << "MONOREPAIR " << goofy << std::endl;
            if (column == m)
                return RepairResult();
            interval_t muUI;
            interval_t muLI;
            interval_t mlUI;
            interval_t mlLI;
            std::pair<RepairInput, RepairInput> r1 = createResults(muc);
            //~ std::cout << "muc is final? " << muc.isFinal() << std::endl;
            //~ std::cout << "muU: isFinal: " << r1.first.isFinal() << ", trueType: " << r1.first.getContent().trueType << std::endl;
            //~ std::cout << "muL: isFinal: " << r1.second.isFinal() << ", trueType: " << r1.second.getContent().trueType << std::endl;
            
            std::pair<RepairInput, RepairInput> r2 = createResults(mlc);
            RepairInput muU = r1.first;
            RepairInput muL = r1.second;
            RepairInput mlU = r2.first;
            RepairInput mlL = r2.second;
            
            lpb_coefficient_t sum = information.getSumFrom(column + 2);

            if (muc.isFinal()
                    //~ BUGFIX second disjunct necessary due to cutTree option
		|| (muc.getContent().node->getUpperChild() == 0 || muc.getContent().node->getLowerChild() == 0)) 
	      {
                if (muc.getContent().trueType == 0) {
                    muUI = interval_t(sum, coeff_pos_infinity);
                    muLI = interval_t(sum, coeff_pos_infinity);
                } else {
                    muUI = interval_t(coeff_neg_infinity, 0);
                    muLI = interval_t(coeff_neg_infinity, 0);
                }
            } else {
                muUI = information.getInterval(*muc.getContent().node->getUpperChild());
                muLI = information.getInterval(*muc.getContent().node->getLowerChild());
            }
	    //std::cout << "muUI: " << muUI << std::endl;
	    //std::cout << "muLI: " << muLI << std::endl;

            if (mlc.isFinal()
                    //~ BUGFIX second disjunct necessary due to cutTree option
                    || (mlc.getContent().node->getUpperChild() == 0 || mlc.getContent().node->getLowerChild() == 0)) {
                if (mlc.getContent().trueType == 0) {
                    mlUI = interval_t(sum, coeff_pos_infinity);
                    mlLI = interval_t(sum, coeff_pos_infinity);
                } else {
                    mlUI = interval_t(coeff_neg_infinity, 0);
                    mlLI = interval_t(coeff_neg_infinity, 0);
                }
            } else {
                mlUI = information.getInterval(*mlc.getContent().node->getUpperChild());
                mlLI = information.getInterval(*mlc.getContent().node->getLowerChild());
            }
	    //std::cout << "mlUI: " << mlUI << std::endl;
	    //std::cout << "mlLI: " << mlLI << std::endl;

	    //Determine the bounds for sigma 
            RepairInput muPrimeSigma = muU;
            RepairInput mlPrimeSigma = mlU;
	    lpb_coefficient_t tmp;


	    if (goofy) {
	    //Determine the tightest upper bound for sigma
	      lpb_coefficient_t u_sigma = muUI.second - mlUI.first;
	      //std::cout << "Initial u_sigma = " << u_sigma << std::endl;

	      tmp = muLI.second - mlLI.first;
	      if (tmp < u_sigma) {
                u_sigma = tmp;
                muPrimeSigma = muL;
                mlPrimeSigma = mlL;
	      }
	      //std::cout << "Final u_sigma = " << u_sigma << std::endl;
	      //Conflict on sigma
	      if (oldbound >= u_sigma) {
		return monorepair(muPrimeSigma, mlPrimeSigma, column + 1, m, oldbound, goofy);
	      } 
	    } else { //non-goofy 
	      //Determine the tightest lower bound for sigma
	      lpb_coefficient_t l_sigma = muUI.first - mlUI.second;
	      //std::cout << "Initial l_sigma = " << l_sigma << std::endl;

	      tmp = muLI.first - mlLI.second;
	      if (tmp > l_sigma) {
                l_sigma = tmp;
		//std::cout << "l_sigma = " << l_sigma << std::endl;
                muPrimeSigma = muL;
                mlPrimeSigma = mlL;
	      }
	      //std::cout << "Final l_sigma = " << l_sigma << std::endl;
	      //Conflict on sigma
	      if (l_sigma >= oldbound) {
		return monorepair(muPrimeSigma, mlPrimeSigma, column + 1, m, oldbound, goofy);
	      }
	    }
	    	    
            // non-recursive part
	    //cut-down version of repair, maybe one would do this differently from scratch
            lpb_coefficient_t ak2 = information.getCoefficient(column + 1);
	    //std::cout << "old coefficient: " << ak2 << std::endl;
            lpb_coefficient_t mus, mlb, nub, nls;
            lpb_coefficient_t leftA = 0;
            lpb_coefficient_t rightA = 0;

	    if (goofy) {
	      if (muUI.second < (muLI.second + ak2)) {
                nub = muUI.second;
	      }
	      else {
                nub = muLI.second;
                ++rightA;
	      }
            
	      //std::cout << "nub: " << nub << std::endl;
	      //std::cout << "rightA: " << rightA << std::endl;

	      if (mlUI.first > (mlLI.first + ak2)) {
                nls = mlUI.first;
	      }
	      else {
                nls = mlLI.first;
                --rightA;
	      }
	      mus = oldbound;
	      mlb = 0;
	      //std::cout << "nls: " << nls << std::endl;
	      //std::cout << "rightA: " << rightA << std::endl;
	    } else {//non-goofy

	      if (muUI.first > (muLI.first + ak2)) {
                mus = muUI.first;
	      }
	      else {
                mus = muLI.first;
                ++leftA;
	      }
	      //std::cout << "mus: " << mus << std::endl;
	      //std::cout << "leftA: " << leftA << std::endl;

	      if (mlUI.second < (mlLI.second + ak2)) {
                mlb = mlUI.second;
	      }
	      else {
                mlb = mlLI.second;
                --leftA;
	      }
	      //std::cout << "mlb: " << mlb << std::endl;
	      //std::cout << "leftA: " << leftA << std::endl;

	      nub = oldbound;
	      nls = 0;
	    }

            const lpb_coefficient_t lhs = mus - mlb - nub + nls;
            const lpb_coefficient_t rhs = rightA - leftA;
            
	    //std::cout << "lhs: " << lhs << std::endl;
	    //std::cout << "rhs: " << rhs << std::endl;

            // TODO: correct?!
            assert(rhs != 0);
            const lpb_coefficient_t bound = lhs / rhs;
	    if (bound <= 0) 
	      {
		//std::cout << "BOUND <= 0" << std::endl;
	      }
	    // if (rhs < 0) {//is upper bound
	    //   return RepairResult(true, !(rhs < 0), ceil(bound), column + 1); } else {
	    //   return RepairResult(true, !(rhs < 0), floor(bound), column + 1);
	    // }


	    return RepairResult(true, !(rhs < 0), bound, column + 1);
	}


    private:
        ChooserT chooser;
        SolverInformation information;
};


} // NAMESPACE_BFR

#endif  // EXTENDEDTREESOLVER_HPP_
